!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE MOD_ICE2D
#  if defined (ICE)
  USE ALL_VARS
  USE MOD_PREC
#  if defined (SPHERICAL)
  USE MOD_SPHERICAL
  !#  if defined (NORTHPOLE)
  USE MOD_NORTHPOLE
  !#  endif
#  endif
  
# if defined (MULTIPROCESSOR)
  USE MOD_PAR
# endif
  
  !===============================================================================|
  !   use CICE      
  !===============================================================================|
  USE ICE_KINDS_MOD
  USE ICE_STATE
  USE ICE_MECHRED
  USE ICE_CONSTANTS
  USE ICE_CALENDAR, only: dtice,dyn_dt,ndyn_dt,dtei
  
  IMPLICIT NONE
  SAVE
  
  Logical :: EVP_ON   
! afm 20210527
  Logical :: basalstress   
  REAL(SP), PARAMETER   :: ECC = 4.0_SP
  REAL(SP), PARAMETER   :: ECCI = C1I/ECC
  real (kind=dbl_kind) :: &
       Tdte             &! subcycling timestep for EVP dynamics, s
       ,  tdamp2           ! 2(wave damping time scale T)
  
  REAL(kind=dbl_kind)   ::  &
       dte2T      &  ! dte/2T
       ,denom1     &  ! constants for stress equation
       ,  denom2        !
  
  real(SP), ALLOCATABLE :: strocnxE(:),strocnyE(:)! ice-ocean stress
  real(SP), ALLOCATABLE :: STRAIRXE(:),STRAIRYE(:)
! afm 20200505 basal stress  
!Yoyo 20200604 for output taubx tauby
  real(SP), ALLOCATABLE :: Tbu(:),taubx(:),tauby(:),vice_elem(:)
 
  REAL(SP), ALLOCATABLE :: AICETMP(:),AIU(:)
  REAL(SP), ALLOCATABLE :: UICE2(:),VICE2(:),UICE2F(:),VICE2F(:)

  REAL(SP), ALLOCATABLE :: UICE2N(:), VICE2N(:)

  REAL(SP), ALLOCATABLE :: UICE2RK(:),VICE2RK(:)
  REAL(SP), ALLOCATABLE :: PSTXICE(:),PSTYICE(:)
  REAL(SP), ALLOCATABLE :: ADVSTRX(:),ADVSTRY(:)
  REAL(SP), ALLOCATABLE :: ADVUICE(:),ADVVICE(:)
  REAL(SP), ALLOCATABLE :: IMASS1(:),IMASS(:),PICE(:)
  REAL(SP), ALLOCATABLE :: TAUXWI(:),TAUYWI(:),VREL(:)
  INTEGER,  ALLOCATABLE :: ISICEN(:),ISICEC(:),ISICEN_OLD(:),ISICEC_OLD(:)

  INTEGER,  ALLOCATABLE :: NN_STRESS(:),NN_STRESS_OLD(:)




#  if defined (ICE_EMBEDDING)  
!----J. Ge for ice baroclinic pressure gradient-----!
  REAL(SP) :: NET_DRHOX ! net baroclinic gradient force for whole ice
  REAL(SP) :: NET_DRHOY ! net baroclinic gradient force for whole ice
  integer,  allocatable :: iwbnd(:) ! cell boundary between ice and water 
  integer,  allocatable :: Kzice(:)  ! bottom sigma index of icepack
  REAL(SP), ALLOCATABLE :: DRHOX_ICE (:,:) !baroclinic gradient force
  REAL(SP), ALLOCATABLE :: DRHOY_ICE (:,:) !baroclinic gradient force
!----J. Ge for ice baroclinic pressure gradient-----!
#  endif

  real (kind=dbl_kind), dimension (:,:) ,allocatable,save ::&     
!     &   shear       ! strain rate II component (1/s)
       divu       &   ! strain rate I component, velocity divergence (1/s)
       ,Delta         ! function of strain rates (1/s)
  
  real (kind=SP), dimension (:) ,allocatable,save ::&
       DivuN        & ! strain rate I component, velocity divergence (1/s)
       ,DeltaN      & ! function of strain rates (1/s)
       ,DivuE       & ! strain rate I component, velocity divergence (1/s)
       ,DeltaE        ! function of strain rates (1/s)
  
  integer (kind=int_kind) ::  &
        kdyn        & ! type of dynamics ( 1 = evp )
      , ndte          ! number of subcycles:  ndte=dyn_dt/dte
  
  REAL(DP), PARAMETER :: eyc =0.36
  REAL(DP) :: rcon
  
  CHARACTER(LEN=80) :: RHEOLOGY

  REAL(SP), PARAMETER :: DRAGW = 0.0055_SP*RHOW
  ! drag coefficient for water on ice *rhow (kg/m^3)
# if defined (ICE_FRESHWATER)
  ! afm 20151112 & EJA 20160921 - turning angle 0
  REAL(SP), PARAMETER :: COSW = C1I  !cos(ocean turning angle) !turning angle = 0
  REAL(SP), PARAMETER :: SINW = C0I  !sin(ocean turning angle) !turning angle = 0
# else
  REAL(SP), PARAMETER :: COSW = 0.9063_SP  !cos(ocean turning angle) !turning angle = 25.
  REAL(SP), PARAMETER :: SINW = 0.4226_SP  !sin(ocean turning angle) !turning angle = 25.
# endif

# if defined (ICE_FRESHWATER)
  ! afm 20151112 & EJA 20160921 - turning angle 0
  REAL(SP), PARAMETER :: COSA = c1i  !cos(wind turning angle) !turning angle = 0.
  REAL(SP), PARAMETER :: SINA = c0i  !sin(wind turning angle) !turning angle = 0.
# else
  REAL(SP), PARAMETER :: COSA = 0.9063_SP  !cos(wind turning angle) !turning angle = 25.
  REAL(SP), PARAMETER :: SINA = 0.4226_SP  !sin(wind turning angle) !turning angle = 25.
# endif

  ! afm 20151112 & EJA 20160921 - not used for freshwater
# if !defined (ICE_FRESHWATER) 
  REAL(SP), PARAMETER :: COSA30 = 0.86602_SP  !cos(wind turning angle) !turning angle = 30.
  REAL(SP), PARAMETER :: SINA30 = 0.50000_SP  !sin(wind turning angle) !turning angle = 30.
 
  REAL(SP), PARAMETER :: COSA20 = 0.93969_SP  !cos(wind turning angle) !turning angle = 20.
  REAL(SP), PARAMETER :: SINA20 = 0.34202_SP  !sin(wind turning angle) !turning angle = 20.
# endif

  !!---------------------------------------------------------------------------------   
  real (kind=dbl_kind), dimension (:) ,allocatable,save ::&
       prs_sig       ! replacement pressure, for stress calc
  REAL(SP), ALLOCATABLE :: SIG1(:), SIG2(:) 
  REAL(SP), ALLOCATABLE :: SIG11(:),SIG22(:),SIG12(:)
  REAL(SP), ALLOCATABLE :: PSIG1(:), PSIG2(:) !!principal_stress
  
  CONTAINS

!===============================================================================|
!===============================================================================|
   SUBROUTINE ALLOC_UVICE
   
   IMPLICIT NONE
  
   CHARACTER(LEN=120) :: FNAME
   INTEGER  :: ISCAN
   
   ALLOCATE(UICE2(0:NT))     ;    UICE2 = 0.0_SP
   ALLOCATE(VICE2(0:NT))     ;    VICE2 = 0.0_SP
   ALLOCATE(AIU(0:NT))       ;      AIU = 0.0_SP
 
   ALLOCATE(strocnxE(0:NT))  ; strocnxE =0.0_SP
   ALLOCATE(strocnyE(0:NT))  ; strocnyE =0.0_SP

   ALLOCATE(ISICEC(0:NT))    ; ISICEC     = 0
   ALLOCATE(ISICEC_OLD(0:NT)); ISICEC_OLD = 0
   ALLOCATE(ISICEN(0:MT))    ; ISICEN     = 0
   ALLOCATE(ISICEN_OLD(0:MT)); ISICEN_OLD = 0

   ALLOCATE(NN_STRESS(0:MT));     NN_STRESS     = 0
   ALLOCATE(NN_STRESS_OLD(0:MT)); NN_STRESS_OLD = 0

#  if defined (ICE_EMBEDDING)  
!----J. Ge for ice baroclinic pressure gradient-----!
   allocate(iwbnd(0:nt))     ;  iwbnd = 0
   allocate(kzice(0:nt))     ;  kzice  = 0
   allocate(DRHOX_ICE(0:nt,kb)) ; DRHOX_ICE = 0.0_sp
   allocate(DRHOY_ICE(0:nt,kb)) ; DRHOY_ICE = 0.0_sp
!----J. Ge for ice baroclinic pressure gradient-----!
#  endif

   !  for ridge redistribution
   allocate(Divu(M,1))       ;  Divu  =0.0
   allocate(Delta(M,1))      ;  Delta =0.0

   allocate(DivuN(0:MT))     ; DivuN  =0.0
   allocate(DeltaN(0:MT))    ; DeltaN =0.0

   allocate(DivuE(0:NT))     ; DivuE  =0.0
   allocate(DeltaE(0:NT))    ; DeltaE =0.0

   ALLOCATE(SIG1(0:MT))      ;  SIG1 =0.0_SP
   ALLOCATE(SIG2(0:MT))      ;  SIG2 =0.0_SP
   ALLOCATE(SIG11(0:MT))     ;  SIG11=0.0_SP
   ALLOCATE(SIG22(0:MT))     ;  SIG22=0.0_SP
   ALLOCATE(SIG12(0:MT))     ;  SIG12=0.0_SP

   ALLOCATE(PRS_SIG(0:MT))   ;  PRS_SIG=0.0_SP
   ALLOCATE(PSIG1(0:MT))     ;  PSIG1  =0.0_SP
   ALLOCATE(PSIG2(0:MT))     ;  PSIG2  =0.0_SP

      dyn_dt=dtice/real(ndyn_dt,kind=dbl_kind)

# if defined (ICE_FRESHWATER)
! afm 20151113 & EJA 20160921
! ndte=120-->30, based on time stepping analysis, it should be allowable
! and saves some computational time
     ndte = 30
# else
     ndte =120 !!  
# endif
 
      Tdte = dyn_dt/real(ndte,kind=dbl_kind)        ! s
      dtei = c1i/Tdte              ! 1/s
      tdamp2 = c2I*eyc*dyn_dt         ! s

      ! constants for stress equation
      dte2T = Tdte/tdamp2                    ! unitless
      denom1 = c1i/(c1i+dte2T)
      denom2 = c1i/(c1i+dte2T*ecc)

      rcon  = 1230._dbl_kind*eyc*dyn_dt*dtei**2  ! kg/s

      IF(MSR) write(ipt,*)' ICE dynamice time step  dyn_dt =' , dyn_dt
   
   RETURN
   END SUBROUTINE ALLOC_UVICE
!===============================================================================|
!
!===============================================================================|
   SUBROUTINE ICE_RUNGE_KUTTA

   IMPLICIT NONE
!   integer (kind=int_kind), intent(in) :: kstrngth
   integer (kind=int_kind) :: kstrngth
   INTEGER :: K   ,KID
   integer :: I,J,I1
   real, allocatable :: rsub(:,:)
   real(SP) :: aice_tmp
   REAL(SP)  ::  UI_TMP,VI_TMP
   INTEGER   ::  N_TMP
   INTEGER   :: IUVN
   REAL(SP)  :: ART_TMP
!===============================================================================|
   !  temp use for calculate the air -ice stress
   real :: vmag_tmp      &! surface wind magnitude   (m/s)
      ,      rd_tmp      &! sqrt of exchange coefficient (momentum)
      ,     rdn_tmp      &! sqrt of exchange coefficient (momentum)
      ,     tau_tmp      &! stress at zlvl
      ,   ustar_tmp      &! ustar (m/s)
      ,    rhoa_tmp      &! air density (kg/m^3)
      ,windspeed_tmp

!===============================================================================|
   ! major/minor axis length ratio, squared
      
   ALLOCATE(AICETMP(0:MT));   AICETMP = 0.0_SP
   ALLOCATE(UICE2F(0:NT)) ;   UICE2F  = 0.0_SP
   ALLOCATE(VICE2F(0:NT)) ;   VICE2F  = 0.0_SP
   ALLOCATE(UICE2RK(0:NT));   UICE2RK = 0.0_SP
   ALLOCATE(VICE2RK(0:NT));   VICE2RK = 0.0_SP
   ALLOCATE(PSTXICE(0:NT));   PSTXICE = 0.0_SP
   ALLOCATE(PSTYICE(0:NT));   PSTYICE = 0.0_SP
   ALLOCATE(ADVSTRX(0:NT));   ADVSTRX = 0.0_SP
   ALLOCATE(ADVSTRY(0:NT));   ADVSTRY = 0.0_SP
   ALLOCATE(ADVUICE(0:NT));   ADVUICE = 0.0_SP
   ALLOCATE(ADVVICE(0:NT));   ADVVICE = 0.0_SP
   ALLOCATE(IMASS1(0:NT)) ;   IMASS1  = 0.0_SP
   ALLOCATE(PICE(0:MT))   ;   PICE    = 0.0_SP
   ALLOCATE(IMASS(0:MT))  ;   IMASS   = 0.0_SP
   ALLOCATE(TAUXWI(0:NT)) ;   TAUXWI  = 0.0_SP
   ALLOCATE(TAUYWI(0:NT)) ;   TAUYWI  = 0.0_SP
   ALLOCATE(VREL(0:NT))   ;   VREL  = 0.0_SP

   ALLOCATE(UICE2N(0:MT)); UICE2N =0.0_SP
   ALLOCATE(VICE2N(0:MT)); VICE2N =0.0_SP

   !-----------------------------------------------------------------
   ! pressure and forcing terms; set sigma=0 for no ice;
   ! initialize ice velocity in cells previously empty to ocn current
   !-----------------------------------------------------------------

   CALL ICE_STRENGTH(KSTRENGTH)
   DO I=1,M !T
      PICE(I) = STRENGTH(I,1)
   END DO
   PICE =PICE*RAMP
# if defined (MULTIPROCESSOR)
   IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,PICE)
# endif

   !-----------------------------------------------------------------
   ! total mass of ice and snow, centered in T-cell
   ! NOTE: vice and vsno must be up to date in all grid cells,
   !       including ghost cells
   !-----------------------------------------------------------------

    ISICEC_OLD = ISICEC  !  ice in last step
    ISICEC = 0
    ISICEN = 0
    IMASS  =0.0_SP
    IMASS1 =0.0_SP

   DO I=1,M 
      IF (AICE(I,1) >p001) then
      IMASS(I) = (RHOI*VICE(I,1) + RHOS*VSNO(I,1)) ! kg/m^2
      IF(IMASS(I) >P01)  ISICEN(I) =1
      END IF
   END DO
# if defined (MULTIPROCESSOR)
      if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,IMASS)
      IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,IMASS)
# endif

   CALL  N2E2D(IMASS,IMASS1)

   DO I=1,M !T 
     AICETMP(I) = AICE(I,1)
   END DO
   AIU =0.0_SP

# if defined (MULTIPROCESSOR)
       IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,AICETMP)
# endif

   CALL N2E2D(AICETMP,AIU)

     !EXCHANGE ELEMENT-BASED VALUES ACROSS THE INTERFACE
#    if defined (MULTIPROCESSOR)
     IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,IMASS1,AIU)
#    endif
   
   DO I=1,N !T
      IF(sum(ISICEN(NV(I,1:3))) > 1 .and.IMASS1(I)>p01  &
         .and. AIU(I)>P001) &
      ISICEC(i)=1
   END DO

#    if defined (MULTIPROCESSOR)
     IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,ISICEC)
#    endif

   DO I=1,N 
     UICE2(I) =UICE2(I)*ISICEC(I)
     VICE2(I) =VICE2(I)*ISICEC(I)
   END DO
#    if defined (MULTIPROCESSOR)
     IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,UICE2,VICE2)
#    endif

   !-----------------------------------------------------------   
   !  calculate the wind stress on ice
   !-----------------------------------------------------------
     allocate(STRAIRXE(0:NT),STRAIRYE(0:NT))
     STRAIRXE = 0.0_SP
     STRAIRYE = 0.0_SP

      !------------------------------------------------------------
      ! Define functions
      !------------------------------------------------------------
      !------------------------------------------------------------
      ! Compute turbulent flux coefficients, wind stress, and 
      !------------------------------------------------------------
      !   rdn_tmp  = vonkar/log(zref/iceruf) ! neutral coefficient
!        rdn_tmp = 0.4/log(10./0.0005)
         DO I=1,N
          windspeed_tmp=SQRT(uuwind(i)*uuwind(i)+vvwind(i)*vvwind(i))
!           if ( aicen(i,j,ni) > puny) then
      !------------------------------------------------------------
      ! define some needed variables
      !------------------------------------------------------------
         vmag_tmp = max(1.0, windspeed_tmp)

# if defined (ICE_FRESHWATER)
! afm turning angle 0 ------------- EJA 20160921
         STRAIRXE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*COSA - VVWIND(I)*SINA)*AIU(I)
         STRAIRYE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*SINA + VVWIND(I)*COSA)*AIU(I)
# else
         IF(vmag_tmp<15.0_SP) THEN
         STRAIRXE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*COSA30 - VVWIND(I)*SINA30)*AIU(I) 
         STRAIRYE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*SINA30 + VVWIND(I)*COSA30)*AIU(I) 
         ELSE
         STRAIRXE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*COSA20 - VVWIND(I)*SINA20)*AIU(I)
         STRAIRYE(I) =1.3E-3*(1.1+0.04*vmag_tmp)*vmag_tmp*(UUWIND(I)*SINA20 + VVWIND(I)*COSA20)*AIU(I)
         ENDIF
# endif

      END DO

        STRAIRXE=STRAIRXE*RAMP
        STRAIRYE=STRAIRYE*RAMP     


if (.true.) then !dwang
!------------J. Ge for baroclinic pressure gradient force----------------------!
#  if defined (ICE_EMBEDDING)
  
   call find_ice_water_boundary
   call get_kzice
   call baropg_ice

#  endif 
end if !dwang
!------------------------------------------------------------------------------!



!-----------------------------------------------------------
!
!------BEGIN MAIN LOOP OVER EXTERNAL MODEL 4 STAGE RUNGE-KUTTA INTEGRATION-----!
!
      EVP_ON= .TRUE.
! afm 20210527
      basalstress= .TRUE.
!      basalstress= .false.
      
      IF(EVP_ON) THEN

         !CALL ICEOCEAN_STRESS

! afm & YY add basal stress flag 20200506
         allocate(Tbu(0:NT),taubx(0:NT),tauby(0:NT),vice_elem(0:NT))
         Tbu = 0.0_DP; taubx = 0.0_DP; tauby = 0.0_DP; vice_elem = 0.0_DP
!         allocate(Tbu(0:NT), vice_elem(0:NT))
!         Tbu = 0.0_DP; vice_elem = 0.0_DP

         IF (basalstress) THEN
         call basal_stress 
         END IF

# if defined(MULTIPROCESSOR)
            IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,Tbu)
            IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,taubx,tauby) 
# endif
! afm 20200505 basal stress

         DO KID=1,ndte              ! subcycling
            
            UICE2N=0.0_SP;VICE2N=0.0_SP
            DO I=1,M
               IF(ISICEN(I)==1)THEN
                  IUVN=0
                  ART_TMP=0.0_SP
                  DO I1=1, NTVE(I)
                     UICE2N(I) =UICE2N(I)+UICE2(NBVE(I,I1))*ISICEC(NBVE(I,I1))
                     VICE2N(I) =VICE2N(I)+VICE2(NBVE(I,I1))*ISICEC(NBVE(I,I1))
                     IUVN   =IUVN+ISICEC(NBVE(I,I1))
                  END DO
                  IF(IUVN>0) THEN
                     UICE2N(I) = UICE2N(I) /FLOAT(IUVN)
                     VICE2N(I) = VICE2N(I) /FLOAT(IUVN)
                  ENDIF
               ENDIF
            END DO
# if defined(MULTIPROCESSOR)
            IF(PAR)CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,UICE2N,VICE2N)
            IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,UICE2N,VICE2N)
# endif
            
            CALL ICEOCEAN_STRESS
            UICE2RK  = UICE2
            VICE2RK  = VICE2
            CALL ICE_EVP(ADVSTRX,ADVSTRY,KID)
            DO K=1,4
               !CALL ICE_EVP(ADVSTRX,ADVSTRY)
               CALL ICE_UV(K)
               UICE2  = UICE2F
               VICE2  = VICE2F
!                UICE2 = 0.0_SP
!                VICE2 = 0.0_SP

               !EXCHANGE ELEMENT-BASED VALUES ACROSS THE INTERFACE
#     if defined (MULTIPROCESSOR)
               IF(PAR)CALL AEXCHANGE(EC,MYID,NPROCS,UICE2,VICE2)
#     endif
            END DO     !! END RUNGE-KUTTA LOOP
         END DO     !! END ice subcycle LOOP
         
         DO I=1,M
            SIG1(I) =SIG1(I) *ISICEN(I)
            SIG2(I) =SIG2(I) *ISICEN(I)
            SIG12(I) =SIG12(I) *ISICEN(I)
            ISICEN_OLD(I) =ISICEN(I)
         END DO
# if defined(MULTIPROCESSOR)
         IF(PAR)CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,SIG1,SIG2,SIG12)
         IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,SIG1,SIG2,SIG12)
# endif
         
      ELSE
         !==========================================================================
         Divu=0.0_SP;     Delta=0.0_SP
         DivuN=0.0_SP;    DeltaN=0.0_SP
         
         DO I=1,M
            IF(ISICEN(I) ==1) THEN
               DO I1=1, NTVE(I)
                  divuN(I) =divuN(I)+divuE(NBVE(I,I1))
                  DeltaN(I) =DeltaN(I)+DeltaE(NBVE(I,I1))
               END DO
               DivuN(I)  =DivuN(I)/float(NTVE(I))
               DeltaN(I) =DeltaN(I)/float(NTVE(I))
            ENDIF
         END DO
         
      ENDIF !! EVP selection

      Divu(1:M,1) =DivuN(1:M)
      Delta(1:M,1)=DeltaN(1:M)

!==========================================================================
      DEALLOCATE(UICE2N,VICE2N)
      
      DEALLOCATE(AICETMP)
      DEALLOCATE(UICE2F,VICE2F)
      DEALLOCATE(UICE2RK,VICE2RK)
      DEALLOCATE(PSTXICE,PSTYICE)
      DEALLOCATE(ADVSTRX,ADVSTRY)
      DEALLOCATE(ADVUICE,ADVVICE)
      DEALLOCATE(IMASS1,IMASS,PICE)
      DEALLOCATE(TAUXWI,TAUYWI,vrel)
      DEALLOCATE(STRAIRXE,STRAIRYE)
   
!      IF (basalstress) THEN
      DEALLOCATE(Tbu)
      DEALLOCATE(taubx)
      DEALLOCATE(tauby)
      DEALLOCATE(vice_elem)
!      END IF

      RETURN
    END SUBROUTINE ICE_RUNGE_KUTTA
!==========================================================================
!
!==============================================================================|
    SUBROUTINE ICEOCEAN_STRESS
      
      IMPLICIT NONE

      REAL(SP) :: WATERX,WATERY
      INTEGER  :: I
      
      DO I=1,N !T
         IF(ISICEC(I) == 1)THEN
            WATERX = U(I,1)*COSW - V(I,1)*SINW
            WATERY = V(I,1)*COSW + U(I,1)*SINW
            ! (magnitude of relative ocean current)*rhow*drag*aice
            VREL(I) = AIU(i)*DRAGW*SQRT((U(I,1)-UICE2(I))**2+(V(I,1)-VICE2(I))**2)  ! m/s
            ! ice/ocean stress
            TAUXWI(I) = VREL(I)*WATERX ! NOTE this is not the entire
            TAUYWI(I) = VREL(I)*WATERY ! ocn stress term
         END IF
      END DO
      
      RETURN
    END SUBROUTINE ICEOCEAN_STRESS
    !==============================================================================|

! afm basal stress 20200505 
!=======================================================================
! Computes basal stress Tbu coefficients (landfast ice)
!
! Lemieux, J. F., B. Tremblay, F. Dupont, M. Plante, G.C. Smith, D. Dumont
! (2015). 
! A basal stress parameterization form modeling landfast ice, J. Geophys. Res. 
! Oceans, 120, 3157-3173.
!
! Lemieux, J. F., F. Dupont, P. Blain, F. Roy, G.C. Smith, G.M. Flato (2016). 
! Improving the simulation of landfast ice by combining tensile strength and a
! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121.
!
! note: Tbu is a part of the Cb as defined in Lemieux et al. 2015 and 2016.
    !==============================================================================|
    subroutine basal_stress

      IMPLICIT NONE

      INTEGER  :: I

      real (DP) :: &
         hu,  & ! volume per unit area of ice at u location (mean thickness)
         hwu, & ! water depth at u location
         hcu, & ! critical thickness at u location
!         Tbu, & ! coefficient for basal stress (N/m^2)
!         k1 = 8.0_DP , &  ! first free parameter for landfast parametrization 
! afm 20210527 a larger k1 for Great Lakes, based on the preliminary evaluation
         k1 = 16.0_DP , &  ! first free parameter for landfast parametrization 
         k2 = 15.0_DP , &  ! second free parameter (N/m^3) for landfast parametrization 
         alphab = 20.0_DP  ! alphab=Cb factor in Lemieux et al 2015

      ! convert vice values to that at elements 
      call  N2E2D(vice(:,1),vice_elem)

      DO I=1,N !T
         IF(ISICEC(I) == 1)THEN
           
            ! convert quantities to u-location
            hwu = h1(i)
            hu  = vice_elem(i)
            
            ! 1- calculate critical thickness
            ! aiu is aice (on nodes) interpolated onto elements
            hcu = aiu(i) * hwu / k1

            ! 2- calculate basal stress factor                    
            Tbu(i) = k2 * max(0.0_DP,(hu - hcu)) * exp(-alphab * (1.0_DP - aiu(I)))

            ! basal stress for outputs
            taubx(i) = -uice2(i)*Tbu(i) / (sqrt(uice2(i)**2 + vice2(i)**2) + 0.0001_DP)
            tauby(i) = -vice2(i)*Tbu(i) / (sqrt(uice2(i)**2 + vice2(i)**2) + 0.0001_DP )

         END IF
      END DO

      RETURN
    END SUBROUTINE BASAL_STRESS
    !==============================================================================|
    !
    !==============================================================================|
    !     ACCUMLATE FLUXES FOR ICE                                                 |
    !==============================================================================|
    
    SUBROUTINE ICE_UV(K)       
      
      !==============================================================================|

#  if defined (SPHERICAL)
      !#  if defined (NORTHPOLE)
      USE MOD_NORTHPOLE
      !#  endif   
#  endif   

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: K
      REAL(SP), DIMENSION(0:NT) :: RESXICE,RESYICE
      REAL(SP) :: UICE2FT,VICE2FT
      REAL(SP) :: UI,VI
      INTEGER  :: I,J
      real(SP) :: WIND_SP
      real(SP) :: cca,ccb,ab2,cc1,cc2
      real(SP) :: TAUX, TAUY
!afm & YY: for Basal stress Cb 20200506
      real(SP) :: Cb
      real    (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind    ! residual velocity for basal stress (m/s)
      
      !==============================================================================|
      !
      !--ACCUMULATE RESIDUALS FOR ICE EQUATIONS--------------------------------------|
      !
      
      UICE2FT = UICE2F(0)
      VICE2FT = VICE2F(0)

      DO I=1,N
         IF(ISICEC(I) == 1)THEN
            !--UPDATE----------------------------------------------------------------------|
            ! alpha, beta are defined in Hunke and Dukowicz (1997), section 3.2

!afm & YY: for Basal stress Cb 20200506
!            cca = imass1(i)/Tdte + vrel(I) * cosw*ALPHA_RK(K) ! alpha, kg/m^2 s
            Cb = 0.0_SP
            IF (basalstress) THEN
            Cb  = Tbu(I) / (sqrt(UICE2(I)**2 + VICE2(I)**2) + u0) ! for basal stress
!! afm debug
!if (cb > 0.0_SP .and. vice_elem(i)>=0.5_SP ) then
!print*, "cca terms 1:", imass1(i)/Tdte, vrel(i)* cosw*ALPHA_RK(K), cb
!print*, "cca terms 2:",  h1(i), vice_elem(i),alpha_rk(k),k 
!endif
            END IF 

            ! revp = 0 for classic evp, 1 for revised evp
            cca = imass1(i)/Tdte + vrel(I) * cosw*ALPHA_RK(K) + Cb ! alpha, kg/m^2 s
            ccb = (cor(i)*imass1(i) + vrel(I) * sinw)*ALPHA_RK(K)  ! beta,  kg/m^2 s
            ab2 = cca**2 + ccb**2
            
            ! divergence of the internal stress tensor
            !  ADVSTRX(I),ADVSTRY(I)
            
            ! finally, the velocity components

!----------------J. Ge for baroclinic pressure gradient force----------------!
# if !defined (ICE_EMBEDDING)
            cc1 = (ADVSTRX(I) - IMASS1(I)*PSTXICE(I) + TAUXWI(I) +STRAIRXE(I))&
                 &       *ALPHA_RK(K) + imass1(I)/Tdte*UICE2RK(I)
            cc2 = (ADVSTRY(I) - IMASS1(I)*PSTYICE(I) + TAUYWI(I) +STRAIRYE(I))&
                 &        *ALPHA_RK(K)+ imass1(I)/Tdte*VICE2RK(I)

# else
            cc1 = (ADVSTRX(I) - IMASS1(I)*PSTXICE(I) + net_drhox + TAUXWI(I) +STRAIRXE(I))&
                 &       *ALPHA_RK(K) + imass1(I)/Tdte*UICE2RK(I)
            cc2 = (ADVSTRY(I) - IMASS1(I)*PSTYICE(I) + net_drhoy + TAUYWI(I) +STRAIRYE(I))&
                 &        *ALPHA_RK(K)+ imass1(I)/Tdte*VICE2RK(I)
# endif
!----------------J. Ge for baroclinic pressure gradient force----------------! 
           
            UICE2F(I) = (cca*cc1 + ccb*cc2)/ab2              ! m/s
            VICE2F(I) = (cca*cc2 - ccb*cc1)/ab2
            
         END IF  !  end isicec   
      END DO
      
#  if defined (SPHERICAL)
      CALL ICE_UV_XY(K)
#  endif
      
      DO I=1,N
         IF(ISBCE(I) == 1) THEN
            UI= UICE2F(I)*(SIN(ALPHA(I)))**2-VICE2F(I)*SIN(ALPHA(I))*COS(ALPHA(I))
            VI=-UICE2F(I)*SIN(ALPHA(I))*COS(ALPHA(I))+VICE2F(I)*(COS(ALPHA(I)))**2
!            UICE2F(I)=UI
!            VICE2F(I)=VI
         END IF
         UICE2F(I) =UICE2F(I)*ISICEC(I)
         VICE2F(I) =VICE2F(I)*ISICEC(I)
      END DO
      
      VICE2F(0) = VICE2FT
      UICE2F(0) = UICE2FT
      
      RETURN
    END SUBROUTINE ICE_UV
!==============================================================================|
!
!==============================================================================|

#  if defined (SPHERICAL)
    !#  if defined (NORTHPOLE)
    !
    !==============================================================================|
    SUBROUTINE ICE_EVP_XY(PSIG11PX,PSIG12PY,PSIG12PX,PSIG22PY,KI) 
      !------------------------------------------------------------------------------|
      !   USE BCS
      !   USE MOD_OBCS
      
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: KI
      REAL(SP), DIMENSION(0:MT) :: PUPX,PVPX,PUPY,PVPY
      REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
      
      REAL(SP) :: E11,E12,E22
      Real(SP) :: DD_divu,DT_tension,DS_shear
      Real(SP) :: DELT,ZETA,ETA
      Real(SP) :: PDelta
      REAL(SP) :: SIG11IJ,SIG12IJ,SIG22IJ
      
      REAL(DP) :: ELIJ,UIJ,VIJ,EXFLUX
      INTEGER  :: I,J,II,K,I1,IA,IB,J1,J2,JTMP
      
      REAL(DP) :: UIJ_TMP,VIJ_TMP
      REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
      REAL(DP) :: XC_TMP,YC_TMP
      REAL(DP) :: DLTXE_TMP,DLTYE_TMP,DLTXC_TMP,DLTYC_TMP
      REAL(DP) :: STXJ1,STYJ1,STXJ2,STYJ2
      
      REAL(SP),DIMENSION(0:NT) :: PSIG11PX,PSIG12PY
      REAL(SP),DIMENSION(0:NT) :: PSIG12PX,PSIG22PY
      
      !================================================================================|
      !----------INITIALIZE STRESS ARRAY ----------------------------------------------!
      ! ALL OTHER PROCESSORS ESCAPE HERE
      IF (NODE_NORTHPOLE .EQ. 0) RETURN
      
      DO II=1,NPCV
         I = NCEDGE_LST(II)
         IA  = NIEC(I,1)
         IB  = NIEC(I,2)
         IF(IA == NODE_NORTHPOLE)THEN
            PUPX(IA) = 0.0_SP
            PVPX(IA) = 0.0_SP
            PUPY(IA) = 0.0_SP
            PVPY(IA) = 0.0_SP
         END IF
         IF(IB == NODE_NORTHPOLE)THEN  
            PUPX(IB) = 0.0_SP
            PVPX(IB) = 0.0_SP
            PUPY(IB) = 0.0_SP
            PVPY(IB) = 0.0_SP
         END IF
      END DO
      
      !---------ACCUMULATE STRESS BY LOOPING OVER CONTROL VOLUME HALF EDGES----------!
      
      DO II=1,NPCV
         I = NCEDGE_LST(II)
         I1  = NTRG(I)
         IA  = NIEC(I,1)
         IB  = NIEC(I,2)
         
         IF(ISICEN(IA) == 1 .OR. ISICEN(IB) == 1)THEN
            UIJ = UICE2(I1)
            VIJ = VICE2(I1)
            
            IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
               UIJ_TMP = -VIJ*COS(XC(I1)*DEG2RAD)-UIJ*SIN(XC(I1)*DEG2RAD)
               VIJ_TMP = -VIJ*SIN(XC(I1)*DEG2RAD)+UIJ*COS(XC(I1)*DEG2RAD)
               
               VX1_TMP = REARTH * DCOS(YIJE(I,1)*DEG2RAD) * DCOS(XIJE(I,1)*DEG2RAD)&
                    * 2._DP /(1._DP+Dsin(YIJE(I,1)*DEG2RAD))
               VY1_TMP = REARTH * DCOS(YIJE(I,1)*DEG2RAD) * DSIN(XIJE(I,1)*DEG2RAD)&
                    * 2._DP /(1._DP+Dsin(YIJE(I,1)*DEG2RAD))
               
               VX2_TMP = REARTH * DCOS(YIJE(I,2)*DEG2RAD) * DCOS(XIJE(I,2)*DEG2RAD)&
                    * 2._DP /(1._DP+Dsin(YIJE(I,2)*DEG2RAD))
               VY2_TMP = REARTH * DCOS(YIJE(I,2)*DEG2RAD) * DSIN(XIJE(I,2)*DEG2RAD)&
                    * 2._DP /(1._DP+Dsin(YIJE(I,2)*DEG2RAD))
               
               DLTXE_TMP = VX2_TMP-VX1_TMP
               DLTYE_TMP = VY2_TMP-VY1_TMP
               
               EXFLUX = -UIJ_TMP*DLTYE_TMP                  
               IF(IA == NODE_NORTHPOLE) THEN    
                  PUPX(IA) = PUPX(IA)-EXFLUX/ART1(IA)
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  PUPX(IB) = PUPX(IB)+EXFLUX/ART1(IB)
               END IF

               EXFLUX = -VIJ_TMP*DLTYE_TMP                  
               IF(IA == NODE_NORTHPOLE) THEN    
                  PVPX(IA) = PVPX(IA)-EXFLUX/ART1(IA)
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  PVPX(IB) = PVPX(IB)+EXFLUX/ART1(IB)
               END IF
               
               EXFLUX =  UIJ_TMP*DLTXE_TMP  
               IF(IA == NODE_NORTHPOLE) THEN    
                  PUPY(IA) = PUPY(IA)-EXFLUX/ART1(IA)
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  PUPY(IB) = PUPY(IB)+EXFLUX/ART1(IB)
               END IF
               
               EXFLUX =  VIJ_TMP*DLTXE_TMP  
               IF(IA == NODE_NORTHPOLE) THEN    
                  PVPY(IA) = PVPY(IA)-EXFLUX/ART1(IA)
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  PVPY(IB) = PVPY(IB)+EXFLUX/ART1(IB)
               END IF
            END IF
         END IF
      END DO
!!!==================================================================================================
      I = NODE_NORTHPOLE

      IF(.false.) then
         PUPX(I)=0.0_SP
         PUPY(I)=0.0_SP
         PVPX(I)=0.0_SP
         PVPY(I)=0.0_SP
         
         DO J=1,NTVE(I)
            I1=NBVE(I,J)
            JTMP=NBVT(I,J)
            J1=JTMP+1-(JTMP+1)/4*3
            J2=JTMP+2-(JTMP+2)/4*3
            !       X11=0.5_SP*(VX(I)+VX(NV(I1,J1)))
            !       Y11=0.5_SP*(VY(I)+VY(NV(I1,J1)))
            VX1_TMP = REARTH * DCOS(VY(I)*DEG2RAD) * DCOS(VX(I)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(I)*DEG2RAD))
            VY1_TMP = REARTH * DCOS(VY(I)*DEG2RAD) * DSIN(VX(I)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(I)*DEG2RAD))
            
            VX2_TMP = REARTH * DCOS(VY(NV(I1,J1))*DEG2RAD) * DCOS(VX(NV(I1,J1))*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(NV(I1,J1))*DEG2RAD))
            VY2_TMP = REARTH * DCOS(VY(NV(I1,J1))*DEG2RAD) * DSIN(VX(NV(I1,J1))*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(NV(I1,J1))*DEG2RAD))
            
            X11=0.5_SP*(VX1_TMP+VX2_TMP)
            Y11=0.5_SP*(VY1_TMP+VY2_TMP)
            
            !       X22=XC(I1)
            !       Y22=YC(I1)
            
            XC_TMP = REARTH * DCOS(YC(I1)*DEG2RAD) * DCOS(XC(I1)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(YC(I)*DEG2RAD))
            YC_TMP = REARTH * DCOS(YC(I)*DEG2RAD) * DSIN(XC(I1)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(YC(I)*DEG2RAD))
            
            X22=XC_TMP
            Y22=YC_TMP
            
            !       X33=0.5_SP*(VX(I)+VX(NV(I1,J2)))
            !       Y33=0.5_SP*(VY(I)+VY(NV(I1,J2)))
            
            VX1_TMP = REARTH * DCOS(VY(I)*DEG2RAD) * DCOS(VX(I)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(I)*DEG2RAD))
            VY1_TMP = REARTH * DCOS(VY(I)*DEG2RAD) * DSIN(VX(I)*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(I)*DEG2RAD))
            
            VX2_TMP = REARTH * DCOS(VY(NV(I1,J2))*DEG2RAD) * DCOS(VX(NV(I1,J2))*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(NV(I1,J2))*DEG2RAD))
            VY2_TMP = REARTH * DCOS(VY(NV(I1,J2))*DEG2RAD) * DSIN(VX(NV(I1,J2))*DEG2RAD)&
                 * 2._DP /(1._DP+Dsin(VY(NV(I1,J2))*DEG2RAD))

            X33=0.5_SP*(VX1_TMP+VX2_TMP)
            Y33=0.5_SP*(VY1_TMP+VY2_TMP)
            
            UIJ = UICE2(I1)
            VIJ = VICE2(I1)
            UIJ_TMP = -VIJ*COS(XC(I1)*DEG2RAD)-UIJ*SIN(XC(I1)*DEG2RAD)
            VIJ_TMP = -VIJ*SIN(XC(I1)*DEG2RAD)+UIJ*COS(XC(I1)*DEG2RAD)
            
            PUPX(I)=PUPX(I)+UIJ_TMP*(Y11-Y33)
            PUPY(I)=PUPY(I)+UIJ_TMP*(X33-X11)
            PVPX(I)=PVPX(I)+VIJ_TMP*(Y11-Y33)
            PVPY(I)=PVPY(I)+VIJ_TMP*(X33-X11)
            
         END DO
         
         PUPX(I) = PUPX(I)/ART1(I)
         PVPX(I) = PVPX(I)/ART1(I)
         PUPY(I) = PUPY(I)/ART1(I)
         PVPY(I) = PVPY(I)/ART1(I)
      END IF !! 
!!!==================================================================================================
      
      IF(ISICEN(I) == 0)RETURN
      
      E11 = PUPX(I)
      E12 = 0.5_SP*(PUPY(I)+PVPX(I))
      E22 = PVPY(I)
      
      ! divergence  =  e_11 + e_22
      DD_divu = E11+E22
      ! tension strain rate  =  e_11 - e_22
      DT_tension = E11-E22
      ! shearing strain rate  =  e_12
      DS_shear = 2.0_SP*E12  
      
      DELT = DD_divu*DD_divu+ECCI*(DT_tension*DT_tension+DS_shear*DS_shear)
      DELT = SQRT(DELT)
      !  DELT= min(max(DELT,2.E-9),PICE(I)/(2.0*4.E+8))
      !  DELT= min(max(DELT,2.E-8),PICE(I)/(2.0*4.E+8))
      !  For 1D case
      !--------------------------------------------------  
      !  calculate stress for ridging
      
      DELT= max(DELT,2.E-9)
      IF (KI == ndte) THEN 
         DeltaN (I)= DELT 
         divuN  (I)= DD_divu 
      ENDIF
      !DELT= min(max(DELT,2.E-8),PICE(I)/(2.0*4.E+8))
      !--------------------------------------------------  
      !       Pressure/Delta
      PDelta = PICE(I)*1.0E-6/DELT
      !Pdelta =min(Pice(I)*1.E-6/Delt,Rcon*Art1(I)*1.0E-6)
      !!=========================================================================================| 
      !     first step stress equations
      !     computing sigma1 and sigma2   
      !     SIG1 (I) = ( SIG1 (I) + dte2T*(PDelta*DD_divu-Pice(I))*1.0E+6)  * denom1 
      SIG1 (I) = ( SIG1 (I) + dte2T*(PDelta*(DD_divu-DELT))*1.0E+6)*denom1 
      SIG2 (I) = ( SIG2 (I) + dte2T*PDelta*DT_tension      *1.0E+6)*denom2 
      SIG12(I) = ( SIG12(I) + dte2T*Pdelta*0.5_SP*DS_shear *1.0E+6)*denom2
      
      !     recover sigma11 and sigma22
      SIG11(I)=(SIG1(I)+SIG2(I))*0.5_SP
      SIG22(I)=(SIG1(I)-SIG2(I))*0.5_SP
      
      !!=========================================================================================|
      prs_sig(I)=PDelta*DELT*1.E+6
      
      DO II=1,NPE
         I = NPEDGE_LST(II)
         IA = IEC(I,1)
         IB = IEC(I,2)
         
         IF(CELL_NORTHAREA(IA) == 1)THEN
            PSIG11PX(IA) = 0.0_SP
            PSIG12PY(IA) = 0.0_SP
            
            PSIG12PX(IA) = 0.0_SP
            PSIG22PY(IA) = 0.0_SP 
            
            PSTXICE(IA) = 0.0_SP
            PSTYICE(IA) = 0.0_SP
         END IF
         IF(CELL_NORTHAREA(IB) == 1)THEN  
            PSIG11PX(IB) = 0.0_SP
            PSIG12PY(IB) = 0.0_SP
            
            PSIG12PX(IB) = 0.0_SP
            PSIG22PY(IB) = 0.0_SP 
            
            PSTXICE(IB) = 0.0_SP
            PSTYICE(IB) = 0.0_SP
         END IF
      END DO
      
      DO II=1,NPE
         I=NPEDGE_LST(II)
         IA=IEC(I,1)
         IB=IEC(I,2)
         J1=IENODE(I,1)
         J2=IENODE(I,2)
         
         !     IF(ISICEC(IA) == 1 .OR. ISICEC(IB) == 1)THEN
         ELIJ=0.5_SP*(EL(J1)+EL(J2))
         ! no surface tilt contribution ggao 0516
         ELIJ=0.0
         
         SIG11IJ=0.5_SP*(SIG11(J1)+SIG11(J2))
         SIG22IJ=0.5_SP*(SIG22(J1)+SIG22(J2))
         SIG12IJ=0.5_SP*(SIG12(J1)+SIG12(J2))
         
         VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
              * 2._DP /(1._DP+sin(VY(IENODE(I,1))*DEG2RAD))
         VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
              * 2._DP /(1._DP+sin(VY(IENODE(I,1))*DEG2RAD))
         
         VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
              * 2._DP /(1._DP+sin(VY(IENODE(I,2))*DEG2RAD))
         VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
              * 2._DP /(1._DP+sin(VY(IENODE(I,2))*DEG2RAD))

         DLTXC_TMP = VX2_TMP-VX1_TMP
         DLTYC_TMP = VY2_TMP-VY1_TMP
         
         EXFLUX    = SIG11IJ*DLTYC_TMP
         IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
            PSIG11PX(IA) = PSIG11PX(IA) - EXFLUX/ART(IA)
            PSIG11PX(IB) = PSIG11PX(IB) + EXFLUX/ART(IB)
         ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
            PSIG11PX(IA) = PSIG11PX(IA) - EXFLUX/ART(IA)
         ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
            PSIG11PX(IB) = PSIG11PX(IB) + EXFLUX/ART(IB)
         END IF
         
         EXFLUX    = SIG12IJ*DLTXC_TMP
         IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
            PSIG12PY(IA) = PSIG12PY(IA) + EXFLUX/ART(IA)
            PSIG12PY(IB) = PSIG12PY(IB) - EXFLUX/ART(IB)
         ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
            PSIG12PY(IA) = PSIG12PY(IA) + EXFLUX/ART(IA)
         ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
            PSIG12PY(IB) = PSIG12PY(IB) - EXFLUX/ART(IB)
         END IF

         EXFLUX    = SIG12IJ*DLTYC_TMP
         IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
            PSIG12PX(IA) = PSIG12PX(IA) - EXFLUX/ART(IA)
            PSIG12PX(IB) = PSIG12PX(IB) + EXFLUX/ART(IB)
         ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
            PSIG12PX(IA) = PSIG12PX(IA) - EXFLUX/ART(IA)
         ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
            PSIG12PX(IB) = PSIG12PX(IB) + EXFLUX/ART(IB)
         END IF

         EXFLUX    = SIG22IJ*DLTXC_TMP
         IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
            PSIG22PY(IA) = PSIG22PY(IA) + EXFLUX/ART(IA)
            PSIG22PY(IB) = PSIG22PY(IB) - EXFLUX/ART(IB)
         ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
            PSIG22PY(IA) = PSIG22PY(IA) + EXFLUX/ART(IA)
         ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
            PSIG22PY(IB) = PSIG22PY(IB) - EXFLUX/ART(IB)
         END IF
         
         !      ACCUMULATE BAROTROPIC FLUX
         IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
            PSTXICE(IA)=PSTXICE(IA)-GRAV_E(IA)*ELIJ*(DLTYC_TMP/ART(IA))/          &
                 (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
            PSTYICE(IA)=PSTYICE(IA)+GRAV_E(IA)*ELIJ*(DLTXC_TMP/ART(IA))/           &
                 (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
            PSTXICE(IB)=PSTXICE(IB)+GRAV_E(IB)*ELIJ*(DLTYC_TMP/ART(IB))/           &
                 (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
            PSTYICE(IB)=PSTYICE(IB)-GRAV_E(IB)*ELIJ*(DLTXC_TMP/ART(IB))/           &
                 (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
            PSTXICE(IA)=PSTXICE(IA)-GRAV_E(IA)*ELIJ*(DLTYC_TMP/ART(IA))/           &
                 (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
            PSTYICE(IA)=PSTYICE(IA)+GRAV_E(IA)*ELIJ*(DLTXC_TMP/ART(IA))/           &
                 (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
            PSTXICE(IB)=PSTXICE(IB)+GRAV_E(IB)*ELIJ*(DLTYC_TMP/ART(IB))/           &
                 (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
            PSTYICE(IB)=PSTYICE(IB)-GRAV_E(IB)*ELIJ*(DLTXC_TMP/ART(IB))/           &
                 (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         END IF
         !     END IF  !  is ice
      END DO
      
      RETURN
    END SUBROUTINE ICE_EVP_XY
    !==============================================================================|
    !
    !==============================================================================|
    !     ACCUMLATE FLUXES FOR ICE                                                 |
    !==============================================================================|
    
    SUBROUTINE ICE_UV_XY(K)       
      
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: K
      REAL(SP), DIMENSION(0:NT) :: RESXICE,RESYICE
      REAL(SP) :: UICE2FT,VICE2FT
      INTEGER  :: I,II
      
      REAL(SP) :: UICE_TMP,VICE_TMP,STRAIRX_TMP,STRAIRY_TMP
      REAL(SP) :: PSTXICE_TMP,PSTYICE_TMP
      REAL(SP) :: UICE2RK_TMP,VICE2RK_TMP,UICE2F_TMP,VICE2F_TMP
      REAL(SP) :: DLTYC_TMP ,DLTXC_TMP

      real(SP) :: cca,ccb,ab2,cc1,cc2
      real(SP) :: UW_TMP,VW_TMP
      real(SP) :: TAUXWI_tmp,TAUYWI_tmp
!afm & YY: for Basal stress Cb 20200506
      real(SP) :: Cb
      real    (kind=dbl_kind) :: u0 = 5.e-5_dbl_kind    ! residual velocity for basal stress (m/s)
      !==============================================================================|
      !
      !--ACCUMULATE RESIDUALS FOR ICE EQUATIONS--------------------------------------|
      !
      ! ALL OTHER PROCESSORS ESCAPE HERE
      IF (NODE_NORTHPOLE .EQ. 0) RETURN
      
      DO II=1,NP
         I=NP_LST(II)
         IF(ISICEC(I) == 1)THEN
            UICE_TMP = -VICE2(I)*COS(XC(I)*DEG2RAD)-UICE2(I)*SIN(XC(I)*DEG2RAD)
            VICE_TMP = -VICE2(I)*SIN(XC(I)*DEG2RAD)+UICE2(I)*COS(XC(I)*DEG2RAD)
            
            UW_TMP = -V(I,1)*COS(XC(I)*DEG2RAD)-U(I,1)*SIN(XC(I)*DEG2RAD)
            VW_TMP = -V(I,1)*SIN(XC(I)*DEG2RAD)+U(I,1)*COS(XC(I)*DEG2RAD)
            
            TAUXWI_TMP = -TAUYWI(I)*COS(XC(I)*DEG2RAD)-TAUXWI(I)*SIN(XC(I)*DEG2RAD)
            TAUYWI_TMP = -TAUYWI(I)*SIN(XC(I)*DEG2RAD)+TAUXWI(I)*COS(XC(I)*DEG2RAD)
            
            STRAIRX_TMP = -STRAIRYE(I)*COS(XC(I)*DEG2RAD)-STRAIRXE(I)*SIN(XC(I)*DEG2RAD)
            STRAIRY_TMP = -STRAIRYE(I)*SIN(XC(I)*DEG2RAD)+STRAIRXE(I)*COS(XC(I)*DEG2RAD)
            
            PSTXICE_TMP = -PSTYICE(I)*COS(XC(I)*DEG2RAD)-PSTXICE(I)*SIN(XC(I)*DEG2RAD)
            PSTYICE_TMP = -PSTYICE(I)*SIN(XC(I)*DEG2RAD)+PSTXICE(I)*COS(XC(I)*DEG2RAD)
            !
            !--UPDATE----------------------------------------------------------------------|
            !
            ! alpha, beta are defined in Hunke and Dukowicz (1997), section 3.2
            !        cca = imass1(i)/dtice + vrel * cosw         ! alpha, kg/m^2 s
!afm & YY add Cb for basal stress 20200506
            Cb = 0.0_SP
            IF (basalstress) THEN
            Cb  = Tbu(I) / (sqrt(UICE2(I)**2 + VICE2(I)**2) + u0) ! for basal stress
            END IF
!            cca = imass1(i)/Tdte + vrel(I) * cosw*ALPHA_RK(K) ! alpha, kg/m^2 s
            cca = imass1(i)/Tdte + vrel(I) * cosw*ALPHA_RK(K) + Cb ! alpha, kg/m^2 s
            !        ccb = fm(i,j)        + vrel * sinw        ! beta,  kg/m^2 s
            ccb = (cor(i)*imass1(i) + vrel(I) * sinw)*ALPHA_RK(K)  ! beta,  kg/m^2 s
            
            ab2 = cca**2 + ccb**2
            
            ! divergence of the internal stress tensor
            !  ADVSTRX(I),ADVSTRY(I)

            UICE2RK_TMP = -VICE2RK(I)*COS(XC(I)*DEG2RAD)-UICE2RK(I)*SIN(XC(I)*DEG2RAD)
            VICE2RK_TMP = -VICE2RK(I)*SIN(XC(I)*DEG2RAD)+UICE2RK(I)*COS(XC(I)*DEG2RAD)
            
            ! finally, the velocity components
            cc1 = (ADVSTRX(I) - IMASS1(I)*PSTXICE_TMP + TAUXWI_TMP+STRAIRX_TMP)&
                 &       *ALPHA_RK(K) + imass1(I)/Tdte*UICE2RK_TMP
            cc2 = (ADVSTRY(I) - IMASS1(I)*PSTYICE_TMP + TAUYWI_TMP+STRAIRY_TMP)&
                 &        *ALPHA_RK(K)+ imass1(I)/Tdte*VICE2RK_TMP
            
            UICE2F_TMP = (cca*cc1 + ccb*cc2)/ab2              ! m/s
            VICE2F_TMP = (cca*cc2 - ccb*cc1)/ab2
            
            UICE2F(I)  = VICE2F_TMP*COS(XC(I)*DEG2RAD)-UICE2F_TMP*SIN(XC(I)*DEG2RAD)
            VICE2F(I)  = UICE2F_TMP*COS(XC(I)*DEG2RAD)+VICE2F_TMP*SIN(XC(I)*DEG2RAD)
            VICE2F(I)  = -VICE2F(I)

         else 
            UICE2F(I)=0.0_SP
            VICE2F(I)=0.0_SP
         endif
         
      END DO
      
      RETURN
    END SUBROUTINE ICE_UV_XY
    !==============================================================================|
    !==============================================================================|
    !#  endif
    !end if defined (NORTHPOLE)
#  endif
    !end if defined (SPHERICAL)
    
    
    !==============================================================================|
    !   Calculate Advection  for ICE variables                                     |
    !==============================================================================|
    
    SUBROUTINE ICE_ADV_fvcom(var_ice,cuice,cvice)       
      
      !------------------------------------------------------------------------------|
      IMPLICIT NONE
      
      REAL(SP), DIMENSION(0:MT)    :: XFLUX
      REAL(SP), DIMENSION(0:MT)    :: PUPX,PUPY,PVPX,PVPY  
      REAL(DP), DIMENSION(0:MT)    :: PIPX,PIPY
      REAL(SP), DIMENSION(3*(NT))  :: DTIJ 
      REAL(SP), DIMENSION(3*(NT))  :: UVN
      
      REAL(SP), DIMENSION(0:NT)  :: cuice,cvice
      REAL(SP), DIMENSION(0:MT)  :: var_ice,var_ice1
      
      REAL(SP) :: UTMP,VTMP,SITAI,FFD,FF1,X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2,XI,YI
      REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2,UN,TTIME,ZDEP
      REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT,STPOINT1,STPOINT2
      REAL(SP) :: FACT,FM1
      INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
#  if defined (SPHERICAL)
      REAL(DP) :: TY,TXPI,TYPI
      REAL(DP) :: XTMP1,XTMP
      REAL(DP) :: X1_DP,Y1_DP,X2_DP,Y2_DP,XII,YII
      REAL(DP) :: X11_TMP,Y11_TMP,X33_TMP,Y33_TMP
      !# if defined (NORTHPOLE)
      REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
      REAL(DP) :: TXPI_TMP,TYPI_TMP
      !# endif
#  endif
      
      REAL(SP) :: VI1MIN, VI1MAX, VI2MIN, VI2MAX
      REAL(SP) :: XMIN, XMAX
      REAL(SP) :: ART_TMP 
      
      
      !------------------------------------------------------------------------------!
      !
      !--Initialize Fluxes-----------------------------------------------------------!
      !
      XFLUX     = 0.0_SP
      UVN       =0.0_SP
      !
      !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
      !
      DO I=1,NCV
         I1=NTRG(I)
         UVN(I) = cvice(I1)*DLTXE(I) - cuice(I1)*DLTYE(I)
      END DO
      
      !
      !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
      !
      
      PIPX  = 0.0_SP
      PIPY  = 0.0_SP
      DO I=1,M
         DO J=1,NTSN(I)-1
            I1=NBSN(I,J)
            I2=NBSN(I,J+1)
            FFD=0.5_SP*(var_ice(I1)+var_ice(I2))
            FF1=0.5_SP*(var_ice(I1)+var_ice(I2))
            
            PIPX(I)=PIPX(I)+FF1*DLTYTRIE(i,j)
            PIPY(I)=PIPY(I)+FF1*DLTXTRIE(i,j)
            
         END DO
         PIPX(I)=PIPX(I)/ART2(I)
         PIPY(I)=PIPY(I)/ART2(I)
         
      END DO
      
      DO I=1,NCV_I
         
         I1=NTRG(I)
         
         IF(ISICEC(I1)==1) THEN
         UVN(I) = cvice(I1)*DLTXE(I) - cuice(I1)*DLTYE(I)
         
         IA=NIEC(I,1)
         IB=NIEC(I,2)

         FIJ1=var_ice(IA)+DLTXNCVE(I,1)*PIPX(IA)+DLTYNCVE(I,1)*PIPY(IA)
         FIJ2=var_ice(IB)+DLTXNCVE(I,2)*PIPX(IB)+DLTYNCVE(I,2)*PIPY(IB)
         
         VI1MIN=MINVAL(var_ice(NBSN(IA,1:NTSN(IA)-1)))
         VI1MIN=MIN(VI1MIN, var_ice(IA))
         VI1MAX=MAXVAL(var_ice(NBSN(IA,1:NTSN(IA)-1)))
         VI1MAX=MAX(VI1MAX, var_ice(IA))
         VI2MIN=MINVAL(var_ice(NBSN(IB,1:NTSN(IB)-1)))
         VI2MIN=MIN(VI2MIN, var_ice(IB))
         VI2MAX=MAXVAL(var_ice(NBSN(IB,1:NTSN(IB)-1)))
         VI2MAX=MAX(VI2MAX, var_ice(IB))
         IF(FIJ1 < VI1MIN) FIJ1=VI1MIN
         IF(FIJ1 > VI1MAX) FIJ1=VI1MAX
         IF(FIJ2 < VI2MIN) FIJ2=VI2MIN
         IF(FIJ2 > VI2MAX) FIJ2=VI2MAX
         
         UN=UVN(I)
         
         EXFLUX=-UN* &
              ((1.0_SP+SIGN(1.0_SP,UN))*FIJ2+(1.0_SP-SIGN(1.0_SP,UN))*FIJ1)*0.5_SP 

         IF(ISICEN(IA)==0 .and. EXFLUX >0.0_SP) EXFLUX=0.0
         IF(ISICEN(IB)==0 .and. EXFLUX <0.0_SP) EXFLUX=0.0
         
         XFLUX(IA)=XFLUX(IA)+EXFLUX
         XFLUX(IB)=XFLUX(IB)-EXFLUX
         
         END IF  ! just calculate the ice
      END DO
      
#    if defined (SPHERICAL)
      !#    if defined (NORTHPOLE)
      CALL ICE_ADV_XY(XFLUX,PIPX,PIPY,var_ice,cuice,cvice)
      !#    endif
#    endif  
      
# if defined (MULTIPROCESSOR)
      IF(PAR)CALL NODE_MATCH(0,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,XFLUX)
# endif
      
      !--Update ice variables----------------------------------------------------------!
      !
     
    var_ICE1=0.0_SP
# if defined (ICE_FRESHWATER)    
    DO I=1,M
! afm 20151207 & EJA 20160921 - Fix bug.
! 1. ART_TMP: Underestimates a control volume (It only counts the ice-covered
!             cell), resulting in overestimating flux.
! 2. "IF(ART_TMP>0.0)&...": Ice only advects to a cell where ice already exists.
! Put them into a classical advection form: var_ICE(I)=var_ice(I)-XFLUX(I)*dyn_dt/ART1(i)
!
!!       IF(ISICEN(I)==1)THEN
!       ART_TMP=0.0
!       DO J=1,NTVE(I)
!!          ART_TMP=ART_TMP+ISICEC(NBVE(I,J))*ART(NBVE(I,J))
!          ART_TMP=ART_TMP+real(ISICEC(NBVE(I,J)),SP)*ART(NBVE(I,J))
!       END DO
!        IF(ART_TMP>0.0)&
!          var_ICE1(I)=var_ice(I)-XFLUX(I)*dyn_dt/(ART_TMP/3.0_SP)
!!       ENDIF
!       var_ICE(I)=var_ICE1(I)
!       if (VAR_ICE(I) < 0.0)VAR_ICE(I) = c0i
        var_ICE(I)=var_ice(I)-XFLUX(I)*dyn_dt/ART1(i)
    END DO
# else
    DO I=1,M
       IF(ISICEN(I)==1)THEN
       ART_TMP=0.0
       DO J=1,NTVE(I)
          ART_TMP=ART_TMP+ISICEC(NBVE(I,J))*ART(NBVE(I,J))
       END DO
        IF(ART_TMP>0.0)&
          var_ICE1(I)=var_ice(I)-XFLUX(I)*dyn_dt/(ART_TMP/3.0_SP)
       ENDIF
       var_ICE(I)=var_ICE1(I)
       if (VAR_ICE(I) < 0.0)VAR_ICE(I) = c0i
    END DO
# endif      
      
! - begin afm 20171113 for ice open boundary -------------
      XFLUX_ICE = XFLUX
! - end afm 20171113 for ice open boundary -------------
      RETURN
    END SUBROUTINE ICE_ADV_fvcom
    !==============================================================================|
    
    !#  if defined (NORTHPOLE)
#  if defined (SPHERICAL)
    
    !==============================================================================|
    SUBROUTINE ICE_ADV_XY(XFLUX,PIPX,PIPY,var_ice,cuice,cvice)               
      !------------------------------------------------------------------------------|
      
      IMPLICIT NONE
      REAL(SP), DIMENSION(0:MT)     :: XFLUX,XFLUX_ADV
      REAL(DP), DIMENSION(0:MT)     :: PIPX,PIPY
      REAL(SP), DIMENSION(3*(NT))   :: DTIJ 
      REAL(SP) :: XI,YI
      REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2 
      REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF   
      REAL(SP) :: FACT,FM1
      INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
      REAL(SP) :: TXPI,TYPI
      
      REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
      REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
      REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
      REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
      REAL(SP) :: PIPX_TMP,PIPY_TMP
      REAL(SP) :: U_TMP,V_TMP
      REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
      integer :: KK
      
      REAL(SP), DIMENSION(0:NT)  :: cuice,cvice
      REAL(SP), DIMENSION(0:MT)  :: var_ice,var_ice1
      !!  ggao edge calculation
      REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
      REAL(SP) :: VI1MIN, VI1MAX, VI2MIN, VI2MAX
      
      !------------------------------------------------------------------------------!
      !
      !--Initialize Fluxes-----------------------------------------------------------!
      !
      ! ALL OTHER PROCESSORS ESCAPE HERE
      IF (NODE_NORTHPOLE .EQ. 0) RETURN
      
      DO II=1,NPCV
         I = NCEDGE_LST(II)
         IA = NIEC(I,1)
         IB = NIEC(I,2)
         IF(IA == NODE_NORTHPOLE)THEN
            XFLUX(IA) = 0.0_SP
            XFLUX_ADV(IA) = 0.0_SP
         ELSE IF(IB == NODE_NORTHPOLE)THEN  
            XFLUX(IB) = 0.0_SP
            XFLUX_ADV(IB) = 0.0_SP
         END IF
      END DO
      
      !
      !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
      !
      DO II=1,NPCV
         I = NCEDGE_LST(II)
         I1=NTRG(I)
      END DO

      !
      !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
      !
      I = NODE_NORTHPOLE
      
      IF(I==0)  RETURN
      
      DO II=1,NPCV
         I = NCEDGE_LST(II)
         I1=NTRG(I)
         IA=NIEC(I,1)
         IB=NIEC(I,2)
         
         IF((IA <= M .AND. IB <= M) .AND. I1 <= N)THEN

            !       XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
            !       YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
            !!  ggao edge calculation
            XIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
            YIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
            
            XIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
            YIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
            XI_TMP =0.5_SP*(XIJE1_TMP+XIJE2_TMP)
            YI_TMP =0.5_SP*(YIJE1_TMP+YIJE2_TMP)
            
            !!  change end 0220/2008

            IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
               VXA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * COS(VX(IA)*DEG2RAD) &
                    * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
               VYA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * SIN(VX(IA)*DEG2RAD) &
                    * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
               
               VXB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * COS(VX(IB)*DEG2RAD) &
                    * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
               VYB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * SIN(VX(IB)*DEG2RAD) &
                    * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
               
               !         IF(IA == NODE_NORTHPOLE)THEN
               DXA=XI_TMP-VXA_TMP
               DYA=YI_TMP-VYA_TMP
               !         ELSE IF(IB == NODE_NORTHPOLE)THEN
               DXB=XI_TMP-VXB_TMP
               DYB=YI_TMP-VYB_TMP
               !	 END IF
               !       END IF
               
               IF(IA == NODE_NORTHPOLE)THEN
                  PIPX_TMP=-PIPY(IB)*COS(VX(IB)*DEG2RAD)-PIPX(IB)*SIN(VX(IB)*DEG2RAD)
                  PIPY_TMP=-PIPY(IB)*SIN(VX(IB)*DEG2RAD)+PIPX(IB)*COS(VX(IB)*DEG2RAD)       
                  
                  FIJ1=var_ice(IA)+DXA*PIPX(IA)+DYA*PIPY(IA)
                  FIJ2=var_ice(IB)+DXB*PIPX_TMP+DYB*PIPY_TMP
                  
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  PIPX_TMP=-PIPY(IA)*COS(VX(IA)*DEG2RAD)-PIPX(IA)*SIN(VX(IA)*DEG2RAD)
                  PIPY_TMP=-PIPY(IA)*SIN(VX(IA)*DEG2RAD)+PIPX(IA)*COS(VX(IA)*DEG2RAD)
                  
                  FIJ1=var_ice(IA)+DXA*PIPX_TMP+DYA*PIPY_TMP
                  FIJ2=var_ice(IB)+DXB*PIPX(IB)+DYB*PIPY(IB)
               END IF
               
               VI1MIN=MINVAL(var_ice(NBSN(IA,1:NTSN(IA)-1)))
               VI1MIN=MIN(VI1MIN, var_ice(IA))
               VI1MAX=MAXVAL(var_ice(NBSN(IA,1:NTSN(IA)-1)))
               VI1MAX=MAX(VI1MAX, var_ice(IA))
               VI2MIN=MINVAL(var_ice(NBSN(IB,1:NTSN(IB)-1)))
               VI2MIN=MIN(VI2MIN, var_ice(IB))
               VI2MAX=MAXVAL(var_ice(NBSN(IB,1:NTSN(IB)-1)))
               VI2MAX=MAX(VI2MAX, var_ice(IB))
               IF(FIJ1 < VI1MIN) FIJ1=VI1MIN
               IF(FIJ1 > VI1MAX) FIJ1=VI1MAX
               IF(FIJ2 < VI2MIN) FIJ2=VI2MIN
               IF(FIJ2 > VI2MAX) FIJ2=VI2MAX
               
               !      IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
               UIJ_TMP = -cvice(I1)*COS(XC(I1)*DEG2RAD)-cuice(I1)*SIN(XC(I1)*DEG2RAD)
               VIJ_TMP = -cvice(I1)*SIN(XC(I1)*DEG2RAD)+cuice(I1)*COS(XC(I1)*DEG2RAD)
               
               VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)
               VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)
               
               VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)
               VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)
               
               DLTXE_TMP = VX2_TMP-VX1_TMP
               DLTYE_TMP = VY2_TMP-VY1_TMP
               
               
               UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
               EXFLUX_TMP = -UVN_TMP*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                    (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
               
               IF(IA == NODE_NORTHPOLE)THEN
                  XFLUX(IA)=XFLUX(IA)+EXFLUX_TMP 
               ELSE IF(IB == NODE_NORTHPOLE)THEN
                  XFLUX(IB)=XFLUX(IB)-EXFLUX_TMP 
               END IF
            END IF
         END IF
      END DO
# if defined (MULTIPROCESSOR)
      !IF(PAR)CALL NODE_MATCH(0,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,XFLUX)
# endif
            
      RETURN
    END SUBROUTINE ICE_ADV_XY
    !==============================================================================|
    !#  endif
#  endif
    !end if defined (NORTHPOLE)
    
    !==============================================================================|
    SUBROUTINE ICE_EVP(XFLUX,YFLUX,KI)     
      !------------------------------------------------------------------------------|
      !   USE BCS
      !   USE MOD_OBCS
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: KI
      REAL(SP), ALLOCATABLE :: PUPX(:),PVPX(:),PUPY(:),PVPY(:)
      
      REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
      REAL(SP) :: ELIJ,UIJ,VIJ,EXFLUX
      INTEGER  :: I,J,K,I1,IA,IB,J1,J2,JTMP
   
      REAL(SP) :: E11,E12,E22
      Real(SP) :: DD_divu,DT_tension,DS_shear
      Real(SP) :: DELT,ZETA,ETA
      Real(SP) :: PDelta
      REAL(SP) :: SIG11IJ,SIG12IJ,SIG22IJ
      REAL(SP) :: XFLUX(0:NT),YFLUX(0:NT)

      !   REAL(SP), ALLOCATABLE :: UICE2N(:), VICE2N(:)
      REAL(SP) :: XTMP,XTMP1
      INTEGER, ALLOCATABLE :: IUVN(:)
      
      REAL(DP) :: DLTXE_TMP,DLTYE_TMP,DLTXC_TMP,DLTYC_TMP
      
      REAL(SP), DIMENSION(0:NT) :: PSIG11PX,PSIG12PY
      REAL(SP), DIMENSION(0:NT) :: PSIG12PX,PSIG22PY
      
      REAL(SP)  :: SIG_TMP
      INTEGER   :: N_TMP
      
      REAL(SP) :: UI,VI
      !==============================================================================|
      ! First calculate the velocity gradient on the nodes
      !----------INITIALIZE FLUX ARRAY ----------------------------------------------!
      
      ALLOCATE(PUPX(0:MT))     ;PUPX= 0.0_SP
      ALLOCATE(PVPX(0:MT))     ;PVPX= 0.0_SP
      ALLOCATE(PUPY(0:MT))     ;PUPY= 0.0_SP
      ALLOCATE(PVPY(0:MT))     ;PVPY= 0.0_SP
      
      Delta= 0.0_SP
      divu = 0.0_SP
      
!!!!================================================================================================
      DO I=1,NCV_I
         I1  = NTRG(I)
         IA  = NIEC(I,1)
         IB  = NIEC(I,2)
         
         IF(ISICEN(IA) == 1 .OR. ISICEN(IB) == 1)THEN
            IF(ISICEC(I1)==1) then
               UIJ = UICE2(I1)
               VIJ = VICE2(I1)
            ELSE
               UIJ = UICE2N(IA)*ISICEN(IA)+UICE2N(IB)*ISICEN(IB)
               VIJ = VICE2N(IA)*ISICEN(IA)+VICE2N(IB)*ISICEN(IB)
            ENDIF
            
            !         UIJ = UICE2(I1)
            !         VIJ = VICE2(I1)
            
            EXFLUX = -UIJ*DLTYE(I)
            PUPX(IA) = PUPX(IA)-EXFLUX
            PUPX(IB) = PUPX(IB)+EXFLUX
            
            EXFLUX = -VIJ*DLTYE(I)
            PVPX(IA) = PVPX(IA)-EXFLUX
            PVPX(IB) = PVPX(IB)+EXFLUX

            EXFLUX =  UIJ*DLTXE(I)
            PUPY(IA) = PUPY(IA)-EXFLUX
            PUPY(IB) = PUPY(IB)+EXFLUX
            
            EXFLUX =  VIJ*DLTXE(I)
            PVPY(IA) = PVPY(IA)-EXFLUX
            PVPY(IB) = PVPY(IB)+EXFLUX
            !!  slip boundary
            !! for the solid boundary edge
            !!     isonb(i)=1:  node on the solid boundary  
            !!   normal velocity == 0

            if(.false.) then
            IF(ISONB(IA)==1 .and. ISONB(IB) == 1) THEN
               UI= UICE2(I1)*(SIN(ALPHA(I1)))**2-VICE2(I1)*SIN(ALPHA(I1))*COS(ALPHA(I1))
               VI=-UICE2(I1)*SIN(ALPHA(I1))*COS(ALPHA(I1))+VICE2(I1)*(COS(ALPHA(I1)))**2
               
#           if defined (SPHERICAL)
               XTMP1 = VX(IB)-VX(IA)
               XTMP = XTMP1*TPI
               IF(XTMP1 > 180.0)THEN
                  XTMP = -360.0*TPI+XTMP
               ELSE IF(XTMP1 < -180.0)THEN
                  XTMP =  360.0*TPI+XTMP
               END IF
               DLTXE_TMP = XTMP*COS(DEG2RAD*YC(I1))
               DLTYE_TMP = (VY(IB)-VY(IA))*TPI
#          else
               DLTXE_TMP = VX(IB)-VX(IA)
               DLTYE_TMP = VY(IB)-VY(IA)
#          endif
               
               EXFLUX =  -UI*DLTYE_TMP*0.5_SP
               PUPX(IA) = PUPX(IA)-EXFLUX
               PUPX(IB) = PUPX(IB)+EXFLUX
               
               EXFLUX = -VI*DLTYE_TMP*0.5_SP
               PVPX(IA) = PVPX(IA)-EXFLUX
               PVPX(IB) = PVPX(IB)+EXFLUX
               
               EXFLUX =  UI*DLTXE_TMP*0.5_SP
               PUPY(IA) = PUPY(IA)-EXFLUX
               PUPY(IB) = PUPY(IB)+EXFLUX
               
               EXFLUX =  VI*DLTXE_TMP*0.5_SP
               PVPY(IA) = PVPY(IA)-EXFLUX
               PVPY(IB) = PVPY(IB)+EXFLUX
            END IF
            !!   end solid boundary---non-slip
           endif  !! slip solid boundary condition
            
         END IF
      END DO
!!!-----------------------------------------------------------------------------
# if defined(MULTIPROCESSOR)
      IF(PAR)CALL NODE_MATCH(0,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,PUPX,PVPX)
      IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,PUPX,PVPX)
      IF(PAR)CALL NODE_MATCH(0,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,PUPY,PVPY)
      IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,PUPY,PVPY)
# endif

      PUPX = PUPX/ART1
      PVPX = PVPX/ART1
      PUPY = PUPY/ART1
      PVPY = PVPY/ART1
      
      DO I=1,M
         IF(ISICEN(I) == 1)THEN
            E11 = PUPX(I)
            E12 = 0.5_SP*(PUPY(I)+PVPX(I))
            E22 = PVPY(I)
            
            ! divergence  =  e_11 + e_22
            DD_divu = E11+E22
            ! tension strain rate  =  e_11 - e_22
            DT_tension = E11-E22
            ! shearing strain rate  =  e_12
            DS_shear = 2.0_SP*E12  
            
            DELT = DD_divu*DD_divu+ECCI*(DT_tension*DT_tension+DS_shear*DS_shear)
            DELT = SQRT(DELT)
            !
            !--------------------------------------------------  
            !       calculate stress for ridging
            DELT= max(DELT,2.E-9)
            IF (KI == ndte) THEN 
               DeltaN (I)= DELT 
               DivuN  (I)= DD_divu 
            ENDIF

            !--------------------------------------------------  
            PDelta = PICE(I)*1.0E-6/DELT
            
            !!====================================================================================|
            !     first step stress equations
            !     computing sigma1 and sigma2
            !     SIG1 (I) = ( SIG1 (I) + dte2T*(PDelta*DD_divu-Pice(I))*1.0E+6)  * denom1
            SIG1 (I) = ( SIG1 (I) + dte2T*(PDelta*(DD_divu-DELT))*1.0E+6)*denom1
            SIG2 (I) = ( SIG2 (I) + dte2T*PDelta*DT_tension      *1.0E+6)*denom2
            SIG12(I) = ( SIG12(I) + dte2T*Pdelta*0.5_SP*DS_shear *1.0E+6)*denom2
            
            !     recover sigma11 and sigma22
            SIG11(I)=(SIG1(I)+SIG2(I))*0.5_SP
            SIG22(I)=(SIG1(I)-SIG2(I))*0.5_SP
            
            !!===================================================================================|
            prs_sig(I)=PDelta*DELT*1.E+6
            
         ENDIF       
      END DO
      
# if defined(MULTIPROCESSOR)
      IF(PAR)CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,SIG1,SIG2)
      IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,SIG1,SIG2)
      IF(PAR)CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,SIG11,SIG12,SIG22)
      IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,SIG11,SIG12,SIG22)
# endif
      
      PSIG11PX = 0.0_SP
      PSIG12PY = 0.0_SP
      PSIG12PX = 0.0_SP
      PSIG22PY = 0.0_SP 
      
      DO I=1,NE
         IA=IEC(I,1)
         IB=IEC(I,2)
         J1=IENODE(I,1)
         J2=IENODE(I,2)
         
         IF(ISICEC(IA) == 1 .OR. ISICEC(IB) == 1)THEN
            ELIJ=0.5_SP*(EL(J1)+EL(J2))
            ! no surface tilt contribution ggao 0516 
            ELIJ=0.0    
            
            SIG11IJ=0.5_SP*(SIG11(J1)+SIG11(J2))
            SIG22IJ=0.5_SP*(SIG22(J1)+SIG22(J2))
            SIG12IJ=0.5_SP*(SIG12(J1)+SIG12(J2))
            
#      if defined (SPHERICAL)
            !for spherical coordinator and domain across 360^o longitude         
            XTMP  = VX(J2)*TPI-VX(J1)*TPI
            XTMP1 = VX(J2)-VX(J1)
            IF(XTMP1 >  180.0_SP)THEN
               XTMP = -360.0_SP*TPI+XTMP
            ELSE IF(XTMP1 < -180.0_SP)THEN
               XTMP =  360.0_SP*TPI+XTMP
            END IF
            !-------------------------------------------------------------
            EXFLUX    = SIG11IJ*DLTYC(I)
            PSIG11PX(IA) = PSIG11PX(IA) - EXFLUX/ART(IA)
            PSIG11PX(IB) = PSIG11PX(IB) + EXFLUX/ART(IB)
            
            EXFLUX    = SIG12IJ*XTMP*COS(DEG2RAD*YC(IA))
            PSIG12PY(IA) = PSIG12PY(IA) + EXFLUX/ART(IA)
            EXFLUX    = SIG12IJ*XTMP*COS(DEG2RAD*YC(IB))
            PSIG12PY(IB) = PSIG12PY(IB) - EXFLUX/ART(IB)
            
            EXFLUX    = SIG12IJ*DLTYC(I)
            PSIG12PX(IA) = PSIG12PX(IA) - EXFLUX/ART(IA)
            PSIG12PX(IB) = PSIG12PX(IB) + EXFLUX/ART(IB)
            
            EXFLUX    = SIG22IJ*XTMP*COS(DEG2RAD*YC(IA))
            PSIG22PY(IA) = PSIG22PY(IA) + EXFLUX/ART(IA)
            EXFLUX    = SIG22IJ*XTMP*COS(DEG2RAD*YC(IB))       
            PSIG22PY(IB) = PSIG22PY(IB) - EXFLUX/ART(IB)
            
            !      ACCUMULATE BAROTROPIC FLUX
            PSTXICE(IA)=PSTXICE(IA)-GRAV_E(IA)*ELIJ*DLTYC(I)/ART(IA)
            PSTYICE(IA)=PSTYICE(IA)+GRAV_E(IA)*ELIJ*XTMP*COS(DEG2RAD*YC(IA))/ART(IA)
            
            PSTXICE(IB)=PSTXICE(IB)+GRAV_E(IB)*ELIJ*DLTYC(I)/ART(IB)
            PSTYICE(IB)=PSTYICE(IB)-GRAV_E(IB)*ELIJ*XTMP*COS(DEG2RAD*YC(IB))/ART(IB)
            
#      else
            EXFLUX    = SIG11IJ*DLTYC(I)
            PSIG11PX(IA) = PSIG11PX(IA) - EXFLUX/ART(IA)
            PSIG11PX(IB) = PSIG11PX(IB) + EXFLUX/ART(IB)
            
            EXFLUX    = SIG12IJ*DLTXC(I)
            PSIG12PY(IA) = PSIG12PY(IA) + EXFLUX/ART(IA)
            PSIG12PY(IB) = PSIG12PY(IB) - EXFLUX/ART(IB)
            
            EXFLUX    = SIG12IJ*DLTYC(I)
            PSIG12PX(IA) = PSIG12PX(IA) - EXFLUX/ART(IA)
            PSIG12PX(IB) = PSIG12PX(IB) + EXFLUX/ART(IB)
            
            EXFLUX    = SIG22IJ*DLTXC(I)
            PSIG22PY(IA) = PSIG22PY(IA) + EXFLUX/ART(IA)
            PSIG22PY(IB) = PSIG22PY(IB) - EXFLUX/ART(IB)
            
            !      ACCUMULATE BAROTROPIC FLUX
            PSTXICE(IA)=PSTXICE(IA)-GRAV_E(IA)*ELIJ*DLTYC(I)/ART(IA)
            PSTYICE(IA)=PSTYICE(IA)+GRAV_E(IA)*ELIJ*DLTXC(I)/ART(IA)
            
            PSTXICE(IB)=PSTXICE(IB)+GRAV_E(IB)*ELIJ*DLTYC(I)/ART(IB)
            PSTYICE(IB)=PSTYICE(IB)-GRAV_E(IB)*ELIJ*DLTXC(I)/ART(IB)
            
#    endif     
         END IF
      END DO
   
#  if defined (SPHERICAL)
      CALL ICE_EVP_XY(PSIG11PX,PSIG12PY,PSIG12PX,PSIG22PY,KI)
#  endif
      
      XFLUX=0.0_SP
      YFLUX=0.0_SP

      DO I=1,N
         IF(sum(ISICEN(NV(I,1:3)))==3) THEN
            XFLUX(I) =  PSIG11PX(I)+PSIG12PY(I)
            YFLUX(I) =  PSIG12PX(I)+PSIG22PY(I)
         END IF
      END DO
      
      DEALLOCATE(PUPX,PUPY,PVPX,PVPY)
      
      RETURN
    END SUBROUTINE ICE_EVP
    !==============================================================================|
    !=======================================================================
    !  principal_stress - computes principal stress for yield curve
    !
    !
    subroutine principal_stress
      !
      ! Computes principal stresses for comparison with the theoretical
      ! yield curve; northeast values
      !
      integer (kind=int_kind) :: i
      do i=1,M
         !if(prs_sig(i) > puny) then
         if(prs_sig(i) > 1.E-8) then
            Psig1(i)=(p5*(sig1(i)+sqrt(sig2(i)**2+c4I*sig12(i)**2)))/prs_sig(i)
            Psig2(i)=(p5*(sig1(i)-sqrt(sig2(i)**2+c4I*sig12(i)**2)))/prs_sig(i)
         else
            Psig1(i)=1.0E+2
            Psig2(i)=1.0E+2
         endif
      enddo
      
    end subroutine principal_stress
        




#  if defined (ICE_EMBEDDING)
!==============================================================================|
    ! subroutine to find ice-water boundary
    ! yding
    !
    ! J. Ge modified for parallel running and 
    !    ice baroclinic pressure gradient
    !
    subroutine FIND_ICE_WATER_BOUNDARY

    implicit none
    !
    INTEGER :: I,K,IOUTTMP,ierr 
    !
    !

    do i = 1, N
      if(aice(nv(i,1),1)>1.0e-5 .and. aice(nv(i,2),1)<1.0e-5 .and. aice(nv(i,3),1)<1.0e-5) then
        iwbnd(i) = 1
      else if (aice(nv(i,2),1)>1.0e-5 .and. aice(nv(i,1),1)<1.0e-5 .and. aice(nv(i,3),1)<1.0e-5) then
        iwbnd(i) = 1
      else if (aice(nv(i,3),1)>1.0e-5 .and. aice(nv(i,1),1)<1.0e-5 .and. aice(nv(i,2),1)<1.0e-5) then
        iwbnd(i) = 1
      else if (aice(nv(i,1),1)>1.0e-5 .and. aice(nv(i,2),1)>1.0e-5 .and. aice(nv(i,3),1)<1.0e-5) then
        iwbnd(i) = 1
      else if (aice(nv(i,1),1)>1.0e-5 .and. aice(nv(i,3),1)>1.0e-5 .and. aice(nv(i,2),1)<1.0e-5) then
        iwbnd(i) = 1
      else if (aice(nv(i,2),1)>1.0e-5 .and. aice(nv(i,3),1)>1.0e-5 .and. aice(nv(i,1),1)<1.0e-5) then
        iwbnd(i) = 1
      else
        iwbnd(i) = 0
      end if
    end do 
   end subroutine FIND_ICE_WATER_BOUNDARY 





!=============================================================================|
!   calculate the bottom sigma index of embedding icepack
!=============================================================================|
   SUBROUTINE GET_Kzice
   
   IMPLICIT NONE
   INTEGER  :: I,I1,I2,I3,I4,J,K
   REAL(SP) :: DTMP,ice_depth

   if(dbg_set(dbg_sbr)) write(ipt,*) "Start: get_kzice "    

   DO I=1,N

     if(iwbnd(i)==0)cycle

     D1(I) = H1(I)+ELF1(I)
     ice_depth = max(QTHICE(NV(I,1)),QTHICE(NV(I,2)),QTHICE(NV(I,3)))
     Kzice(I) = 0
     DTMP = 0.0_SP

     DO K=1,KB
       DTMP = DTMP + DZ1(I,K)*D1(I)
       IF(DTMP <= ice_depth)THEN
         Kzice(I) = Kzice(I) + 1
       END IF
     END DO  	 
     
   END DO

   if(dbg_set(dbg_sbr)) write(ipt,*) "End: get_kzice " 

   Return

 End SUBROUTINE GET_Kzice


!==============================================================================|
!     CALCULATE THE BAROCLINIC PRESSURE GRADIENT IN SIGMA COORDINATES          |
!==============================================================================|

   SUBROUTINE BAROPG_ICE

!==============================================================================|
   USE ALL_VARS
   USE MOD_SPHERICAL
   USE MOD_NORTHPOLE
   USE MOD_WD

#  if defined (MULTIPROCESSOR)
   USE MOD_PAR
#  endif

   IMPLICIT NONE
   REAL(SP) :: RIJK(0:N,3,KBM1), DRIJK1(0:N,3,KBM1), DRIJK2(0:N,KBM1)
   REAL(SP) :: TEMP,DIJ,DRHO1,DRHO2
   INTEGER  :: I,K,J,J1,J2,IJK
#  if defined (SPHERICAL)
   REAL(SP) :: XTMP,XTMP1
#  endif

   integer :: ierr
   integer,allocatable, dimension(:)   :: iwbndtmp
   real(sp),allocatable,dimension(:)   :: art1tmp
   real(sp),allocatable,dimension(:,:) :: aicetmp,vicetmp,drhoxtmp,drhoytmp
   real(dp)::total_mass,total_drhox,total_drhoy


!==============================================================================|

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: baropg.F"
   
   ! USE RAMP CALCULATED IN 'internal_step.F'

!----------SUBTRACT MEAN DENSITY TO MINIMIZE ROUNDOFF ERROR--------------------!

   RHO1(:,1:KBM1) = RHO1(:,1:KBM1) - RMEAN1(:,1:KBM1)
   RHO = RHO - RMEAN 

!----------INITIALIZE ARRAYS---------------------------------------------------!

   DRHOX_ice  = 0.0_SP
   DRHOY_ice  = 0.0_SP
   RMEAN(0,:) = 0.0_SP
   RHO(0,:)   = 0.0_SP
   RIJK       = 0.0_SP
   DRIJK1     = 0.0_SP
   DRIJK2     = 0.0_SP

!----------CALCULATE AVERAGE DENSITY ON EACH EDGE------------------------------!

   DO K=1,KBM1
     DO I=1,N
!------J. Ge-------
       if(iwbnd(i)==0)cycle
       if(k>kzice(i))cycle
!------J. Ge-------
       DO J=1,3
         J1=J+1-INT((J+1)/4)*3
         J2=J+2-INT((J+2)/4)*3
         RIJK(I,J,K)  = 0.5_SP*(RHO1(NV(I,J1),K)+RHO1(NV(I,J2),K))
       END DO
     END DO
   END DO

   DO I=1,N
!------J. Ge-------
     if(iwbnd(i)==0)cycle
!------J. Ge-------
     DO J=1,3
       DRIJK1(I,J,1)=RIJK(I,J,1)*(-ZZ1(I,1))
       DO K=2,KBM1
!------J. Ge-------
         if(k>kzice(i))cycle
!------J. Ge-------
         DRIJK1(I,J,K)=0.5_SP*(RIJK(I,J,K-1)+RIJK(I,J,K))*(ZZ1(I,K-1)-ZZ1(I,K))
         DRIJK1(I,J,K)=DRIJK1(I,J,K)+DRIJK1(I,J,K-1)
       END DO
     END DO
   END DO

   DO I=1,N
!------J. Ge-------
     if(iwbnd(i)==0)cycle
!------J. Ge-------
     DRIJK2(I,1)=0.0_SP
     DO K=2,KBM1
!------J. Ge-------
         if(k>kzice(i))cycle
!------J. Ge-------
       DRIJK2(I,K)=0.5_SP*(ZZ1(I,K-1)+ZZ1(I,K))*(RHO(I,K)-RHO(I,K-1)) 
       DRIJK2(I,K)=DRIJK2(I,K-1)+DRIJK2(I,K)
     END DO
   END DO

   DO I = 1, N
!------J. Ge-------
     if(iwbnd(i)==0)cycle
!------J. Ge-------
#  if defined (WET_DRY)
    IF(ISWETCT(I)*ISWETC(I) == 1 .AND. &
      (H(NV(I,1)) > STATIC_SSH_ADJ .OR. H(NV(I,2)) > STATIC_SSH_ADJ .OR. H(NV(I,3)) > STATIC_SSH_ADJ))THEN
#  endif
     DO K=1,KBM1
!------J. Ge-------
        if(k>kzice(i))cycle
!------J. Ge-------
        DO J = 1, 3
          J1=J+1-INT((J+1)/4)*3
          J2=J+2-INT((J+2)/4)*3
          IJK=NBE(I,J)
          DIJ=0.5_SP*(DT(NV(I,J1))+DT(NV(I,J2)))

#         if defined (SPHERICAL)
          DRHO1=-DELTUY(I,J)*DRIJK1(I,J,K)*DT1(I)
          DRHO2=-DELTUY(I,J)*DIJ*DRIJK2(I,K)
#         else
          DRHO1=(VY(NV(I,J1))-VY(NV(I,J2)))*DRIJK1(I,J,K)*DT1(I)
          DRHO2=(VY(NV(I,J1))-VY(NV(I,J2)))*DIJ*DRIJK2(I,K)
#         endif
          DRHOX_ice(I,K)=DRHOX_ice(I,K)+DRHO1+DRHO2

#         if defined (SPHERICAL)
          XTMP  = VX(NV(I,J2))*TPI-VX(NV(I,J1))*TPI
          XTMP1 = VX(NV(I,J2))-VX(NV(I,J1))
          IF(XTMP1 >  180.0_SP)THEN
            XTMP = -360.0_SP*TPI+XTMP
          ELSE IF(XTMP1 < -180.0_SP)THEN
            XTMP =  360.0_SP*TPI+XTMP
          END IF  

          DRHO1=XTMP*COS(DEG2RAD*YC(I))*DRIJK1(I,J,K)*DT1(I)
          DRHO2=XTMP*COS(DEG2RAD*YC(I))*DIJ*DRIJK2(I,K)
!          DRHO1=DELTUX(I,J)*DRIJK1(I,J,K)*DT1(I)
!          DRHO2=DELTUX(I,J)*DIJ*DRIJK2(I,K)
#         else
	  DRHO1=(VX(NV(I,J2))-VX(NV(I,J1)))*DRIJK1(I,J,K)*DT1(I)
          DRHO2=(VX(NV(I,J2))-VX(NV(I,J1)))*DIJ*DRIJK2(I,K)
#         endif
          DRHOY_ice(I,K)=DRHOY_ice(I,K)+DRHO1+DRHO2

       END DO
     END DO
#  if defined (WET_DRY)
    END IF
#  endif
   END DO

#  if defined (SPHERICAL)
   CALL BAROPG_XY(DRIJK1,DRIJK2)
#  endif  


!----------MULTIPLY BY GRAVITY AND ELEMENT DEPTH-------------------------------!

   DO K=1,KBM1
     DRHOX_ice(:,K)=DRHOX_ice(:,K)*DT1(:)*DZ1(:,K)*GRAV_E(:)*RAMP
     DRHOY_ice(:,K)=DRHOY_ice(:,K)*DT1(:)*DZ1(:,K)*GRAV_E(:)*RAMP
   END DO

!----------ADD MEAN DENSITY BACK ON--------------------------------------------!

   RHO1 = RHO1 + RMEAN1
   RHO  = RHO  + RMEAN





   IF(SERIAL)THEN
!----------calculate the whole ice mass----------------------------------------!
    total_mass = 0.0
    do i=1,m
      total_mass = total_mass + art1(i)*vice(i,1)*aice(i,1)*rhoi
    end do

!----------accumulate all gradient force---------------------------------------!
    total_drhox = 0.0
    total_drhoy = 0.0
    do i=1,n
       if(iwbnd(i)==0)cycle
       total_drhox = total_drhox + sum(DRHOX_ice (i,:))
       total_drhoy = total_drhoy + sum(DRHOy_ice (i,:))
    end do

   ENDIF


#  if defined (MULTIPROCESSOR) 
   IF(PAR)THEN

     CALL MPI_BARRIER(MPI_FVCOM_GROUP,IERR)

     ALLOCATE(aicetmp(0:mGL,1))  
     allocate(art1tmp(0:mgl)) 
     allocate(vicetmp(0:mgl,1)) 
     allocate(drhoxtmp(0:ngl,kb))
     allocate(drhoytmp(0:ngl,kb))
     allocate(iwbndtmp(0:ngl)) 

     CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,art1,art1tmp)
     CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,vice,vicetmp)
     CALL ACOLLECT(MYID,MSRID,NPROCS,NMAP,aice,aicetmp)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,DRHOX_ice,drhoxtmp)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,DRHOY_ice,drhoytmp)
     CALL ACOLLECT(MYID,MSRID,NPROCS,EMAP,iwbnd,iwbndtmp)

!----------calculate the whole ice mass----------------------------------------!
     total_mass = 0.0
     do i=1,mgl
       total_mass = total_mass + art1tmp(i)*vicetmp(i,1)*aicetmp(i,1)*rhoi
     end do

!----------accumulate all gradient force---------------------------------------!
     total_drhox = 0.0
     total_drhoy = 0.0
     do i=1,ngl
       if(iwbndtmp(i)==0)cycle
       total_drhox = total_drhox + sum(DRHOXtmp(i,:))
       total_drhoy = total_drhoy + sum(DRHOytmp(i,:))
     end do     

   END IF

#  endif

!----------put the total gradient force on whole icepack-----------------------!
   net_DRHOX = total_drhox / total_mass
   net_DRHOY = total_drhoy / total_mass


   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: baropg.F"

   RETURN
   END SUBROUTINE BAROPG_ICE
!==============================================================================|





!==============================================================================|

#  endif


#  endif
  END MODULE MOD_ICE2D
 
